import networkx as nx
import matplotlib.pyplot as plt
import operator
import torch
import scipy as sp
import math

def HITS(G, epsilon):
    N = G.number_of_nodes()
    number_of_iteration = 0
    # Initialize
    a_old = torch.full([N], 1.0 / math.sqrt(N), dtype = torch.float32)
    h_old = a_old
    A = torch.tensor(nx.adjacency_matrix(G).todense(), dtype = torch.float32)
    h_new = torch.mv(A, a_old)
    a_new = torch.mv(A.T, h_old)
    s_h = 0
    s_a = 0
    for i in range(N):
        s_h += h_new[i] ** 2
        s_a += a_new[i] ** 2
    s_h = math.sqrt(s_h)
    s_a = math.sqrt(s_a)
    h_new /= s_h
    a_new /= s_a
        
    # Loops
    while (torch.norm(a_new - a_old) > epsilon) or (torch.norm(h_new - h_old) > epsilon):
        a_old = a_new
        h_old = h_new
        h_new = torch.mv(A, a_old)
        a_new = torch.mv(A.T, h_old)
        # Normalize
        s_h = 0
        s_a = 0
        for i in range(N):
            s_h += h_new[i] ** 2
            s_a += a_new[i] ** 2
        s_h = math.sqrt(s_h)
        s_a = math.sqrt(s_a)
        h_new /= s_h
        a_new /= s_a
        number_of_iteration += 1
    
    return [h_new, a_new, number_of_iteration]